A trivial SOS decomposition example

Contributed by: votroto Adapted from: (Blekherman et al., 2012; Examples 3.25)

using DynamicPolynomials
using SumOfSquares
import CSDP

The polynomial p = x^2 - x*y^2 + y^4 + 1 is SOS. We can, for example, decompose it as p = 3/4*(x - y^2)^2 + 1/4*(x + y)^2 + 1, which clearly proves that p is SOS, and there are infinitely many other ways to decompose p into sums of squares.

We can use SumOfSquares.jl to find such decompositions.

First, setup the polynomial of interest.

@polyvar x y
p = x^2 - x*y^2 + y^4 + 1

\[ 1 + x^{2} - xy^{2} + y^{4} \]

Secondly, constrain the polynomial to be nonnegative. SumOfSquares.jl transparently reinterprets polyonmial nonnegativity as the appropriate SOS certificate for polynomials nonnegative on semialgebraic sets.

model = SOSModel(CSDP.Optimizer)
@constraint(model, cref, p >= 0)

\[ (1) \cdot 1 + (1) \cdot x^{2} + (-1) \cdot xy^{2} + (1) \cdot y^{4} \text{ is SOS} \]

Thirdly, optimize the feasibility problem!

optimize!(model)

CSDP 6.2.0
This is a pure primal feasibility problem.
Iter:  0 Ap: 0.00e+00 Pobj:  0.0000000e+00 Ad: 0.00e+00 Dobj:  0.0000000e+00
Iter:  1 Ap: 9.00e-01 Pobj:  0.0000000e+00 Ad: 1.00e+00 Dobj:  2.7659722e+01
Iter:  2 Ap: 1.00e+00 Pobj:  0.0000000e+00 Ad: 1.00e+00 Dobj:  1.2562719e+01

Lastly, recover a SOS decomposition. In general, SOS decompositions are not unique!

sos_dec = sos_decomposition(cref, 1e-4)

(-0.7423566403018642·1 - 0.5916342145954683·x + 0.9492757372229985·y^2)^2 + (-6.60921131771784e-17·1 + 1.1201589623631873·y - 2.4872525425116695e-16·x - 1.0534655002111817e-16·y^2)^2 + (0.6232480104529257·1 - 0.782024243528594·x + 5.551115123125783e-17·y^2)^2 + (-0.2459035096662894·1 - 1.2342687879752308e-17·y - 0.19597713808894587·x - 0.31444486753600015·y^2)^2

Converting, rounding, and simplifying - Huzza, Back where we began!

polynomial(sos_dec, Float32)

\[ 1.0 - 1.4199712 \cdot 10^{-16}y + 7.4940054 \cdot 10^{-16}x - 5.551115 \cdot 10^{-16}y^{2} - 5.5238587 \cdot 10^{-16}xy + x^{2} - 2.2824759 \cdot 10^{-16}y^{3} - xy^{2} + y^{4} \]

A deeper explanation and the unexplained 1e-4 parameter

p = x^2 - x*y^2 + y^4 + 1 can be represented in terms of its Gram matrix as

gram = gram_matrix(cref)

GramMatrix with row/column basis:
 SubBasis{Monomial}([1, y, x, y^2])
And entries in a 4×4 SymMatrix{Float64}:
  1.0                    0.0                 …  -0.6273780504812863
  0.0                    1.2547561009625725      0.0
  4.933553565678038e-17  0.0                    -0.5
 -0.6273780504812863     0.0                     1.0

gram.basis.monomials' * gram.Q * gram.basis.monomials

\[ 1.0 + 9.867107131356075 \cdot 10^{-17}x + x^{2} - xy^{2} + y^{4} \]

where the matrix gram.Q is positive semidefinite, because p is SOS. If we could only get the decomposition gram.Q = V' * V, the SOS decomposition would simply be ||V * monomials||^2.

Unfortunately, we can not use Cholesky decomposition, since gram Q is only semidefinite, not definite. Hence, SumOfSquares.jl uses SVD decomposition instead and discards small singular values (in our case 1e-4).

This page was generated using Literate.jl.